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The fundamental problem faced in quantum chemistry is the calculation of 
molecular properties, which are of practical importance in fields ranging from 
materials science to biochemistry. Within chemical precision, the total energy 
of a molecule as well as most other properties, can be calculated by solving the 
Schrodinger equation. However, the computational resources required to obtain 
exact solutions on a conventional computer generally increase exponentially with 
the number of atoms involved^ ^. This renders such calculations intractable for 
all but the smallest of systems. Recently, an efficient algorithm has been pro- 
posed enabling a quantum computer to overcome this problem by achieving 
only a polynomial resource scaling with system size ^ ^ ^. Such a tool would 
therefore provide an extremely powerful tool for new science and technology. 
Here we present a photonic implementation for the smallest problem: obtain- 
ing the energies of H2, the hydrogen molecule in a minimal basis. We perform 
a key algorithmic step — the iterative phase estimation algorithmpEEEl — full, 
achieving a high level of precision and robustness to error. We implement other 
algorithmic steps with assistance from a classical computer and explain how this 
non-scalable approach could be avoided. Finally, we provide new theoretical re- 
sults which lay the foundations for the next generation of simulation experiments 
using quantum computers. We have made early experimental progress towards 
the long-term goal of exploiting quantum information to speed up quantum 
chemistry calculations. 

Experimentalists are just beginning to command the level of control over quantum sys- 
tems required to explore their information processing capabilities. An important long-term 
application is to simulate and calculate properties of other many-body quantum systems. Pi- 
oneering experiments were first performed using nuclear-magnetic-resonance-based systems 
to simulate quantum oscillator^, leading up to recent simulations of a pairing Hamilto- 
niarP^. Very recently the phase transitions of a two-spin quantum magnet were simulatecP^ 
using an ion-trap system. Here we simulate a quantum chemical system and calculate its 
energy spectrum, using a photonic system. 

Molecular energies are represented as the eigenvalues of an associated time-independent 
Hamiltonian H and can be efficiently obtained to fixed accuracy, using a quantum algorithm 
with three distinct steps^: encoding a molecular wavefunction into qubits; simulating its time 
evolution using quantum logic gates; and extracting the approximate energy using the phase 
estimation algorithmP^. The latter is a general-purpose quantum algorithm for evaluating 
the eigenvalues of arbitrary Hermitian or unitary operators. The algorithm estimates the 
phase, 0, accumulated by a molecular eigenstate, under the action of the time-evolution 
operator, U=e~^^''^^, i.e., 

where E is the energy eigenvalue of Therefore, estimating the phase for each eigenstate 
amounts to estimating the eigenvalues of the Hamiltonian. 

We take the standard approach to quantum-chemical calculations by solving an approx- 
imate Hamiltonian created by employing the Born-Oppenheimer approximation (where the 
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electronic Hamiltonian is parameterized by the nuclear configuration) and choosing a suit- 
able truncated basis set in which to describe the non-relativistic electronic system. Typical 
sets consist of a finite number of single-electron atomic orbitals, which are combined to form 
antisymmetric multi-electron molecular states (configurations)?^. Calculating the eigenval- 
ues of the electronic Hamiltonian using all configurations gives the exact energy in the basis 
set and is referred to as full configuration interaction (FCI). For orbitals and m electrons 
there are ways to allocate the electrons among the orbitals. This exponential 

growth is the handicap of FCI calculations on classical computers. 

As described in the Methods Summary, the Hamiltonian is block diagonal in our choice of 
basis, with 2x2 sub- matrices [H'^^'^^ and H^^''^''). We map the configurations spanning each 
sub-space to the qubit computational basis. Since the subspaces are two-dimensional, one 
qubit suffices to represent the wavefunction. The corresponding time-evolution operators, 
lj{p,q)—^-tH(P'ih/h — -y^^]2g];g (p, g) = (l,6) or (3,4) — are therefore one-qubit operators. Finding 
the eigenvalues of each separately, using a phase estimation algorithm, amounts to per- 
forming FCI. For the purpose of our demonstration, we encode exact eigenstates, obtained 
via a preliminary calculation on a classical computer. In our Appendix we show that the 
algorithm is in fact robust to imperfect eigenstate encoding. 

We implement the iterative phase estimation algorithm-^- (IPEA), which advantageously 
reduces the number of qubits and quantum logic gates required. Figjl] a shows the IPEA 
at iteration k. The result of a logical measurement of the top 'control' qubit after each 
iteration determines the fc*^ bit of the binary expansion!^ of (p. Let K bits of this expansion 
be 0=O.0i02---0m; such that 0=0+52"^ where 5 is a remainder 0<5<1. An accuracy of 
±2~^ is achieved with error probabilitj^ e<l— S/vr^RiO.lQ, which is independent of K (the 
bound is saturated for 5=0.5). This error can be eliminated by simply repeating each IPEA 
iteration multiple (n) times, yielding n possible values for the corresponding bit; a majority 
vote of these samples determines the result (see Methods, Section C). 

To resolve the energy differences relevant to chemical processed, absolute molecular ener- 
gies must be computed to an accuracy greater than ^10~^i?h {^khT at room temperature). 
Therefore it is important to demonstrate that the IPEA can achieve the necessary phase 
precision of ~16 bits (the accuracy of the non-relativistic Born-Oppenheimer energy is then 
limited only by the choice of basis). We implement the IPEA with a photonic architecture, 
encoding qubits in polarization of single photons, Figjl] c. Our experiment is possible due 
to the recent development of a photonic two-qubit controUed-unitary quantum logic gate, 
which combines multiple photonic degrees of freedom, linear optical elements and projective 
measurement to achieve the required nonlinear interaction between photonJ^. Such gates 
are high-quality, well-characterized, and have several in-principle paths to scalable optical 
quantum computin^^^. We note that our implementation of two-qubit quantum IPEA is 
the first, in any context, to use entangling gates, outside of a liquid-state ensemble NMR 
architecturd^, which is arguably an in-principle non-scalable architecturd^. We remark that 
an implementation of a semiclassical quantum Fourier transform was performed in iond^, 
combining single-qubit measurement and rotations in place of entangling gates. 

Fig. |2] shows our results: H2 energies calculated over a range of internuclear separations, 
thus reconstructing the potential energy surfaces. Each point is obtained using a 20-bit IPEA 
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and employing n=31 samples per bit. In every case, the algorithm successfully returned the 
energy to within the target precision of ±(2~^°x27r)£'h~10~^i?h, in the minimal basis. For 
example, the ground state energy obtained at the equilibrium bond length, 1.3886 oq (where 
oo is the Bohr radius), is —0.20399±0. 00001 -Eh, which agrees exactly with the result obtained 
on a classical computer to an uncertainty in the least significant bit. 

Achieving such high precision will become a far more significant challenge for large-scale 
implementations: due to the small-scale of our demonstration, we are able to implement each 
power of lJ^P''i'> directly, by re-encoding the same number of gates (see Methods Summary). 
Therefore, the probability of error introduced by gate imperfections remains a constant for 
each bit. This is the main algorithmic feature that allows the high precision obtained in this 
experiment. However, following current proposals, the circuit network for JJ will not gener- 
ally have the same form as for larger implementations (detailed in the Appendix, Section 
A). For each additional digit of precision sought, the gate requirements of the algorithm are 
roughly doubled, thereby amplifying any gate error. 

Important next experimental steps are to demonstrate the two parts of the quantum 
algorithm that we implemented with assistance from a classical computer. Firstly, encod- 
ing even low fidelity eigenstate approximations into qubits is a highly non-trivial step for 
molecules much larger than H2. In many cases this problem could be overcome using a 
heuristic adiabatic state preparation algorithmP^^'^I] Here, ground state approximations, 
for example, can be efficiently obtained provided that the energy gap between the ground 
state and excited states is sufficiently large along the path of the adiabatic evolution!^. 
Secondly, directly calculating and decomposing the molecular evolution operator into logic 
gates does not scale efficiently with molecular sizd^ and an alternative scheme must be em- 
ployed. The proposed solution exploits the fact that the general molecular Hamiltonian is 
a sum of fixed-sized one- and two-electron terms that can be individually efficiently simu- 
lated and combined to approximate the global evolutiorP^. We give an overview of this 
'operator-splitting technique' in the Appendix (Section A) and find that the total number 
of elementary quantum gates required to simulate the evolution of an arbitrary molecule, 
without error correction, scales as 0{N^), where N is the number of single-particle basis 
functions used to describe the molecular system. In this scheme, N is also the number of 
qubits necessary. We also present the quantum logic circuits required to simulate each term 
in the general molecular Hamiltonian — these are the building blocks of a universal quantum 
molecular simulator. Finally, we perform an accurate resource count to reproduce our H2 
simulation in this scalable way: 4 qubits and ~522 perfect gates are required to simulate 
the full unitary propagator such that the error of the simulated evolution is within chemical 
precision. 

Other major challenges in the path to scalability include those associated with scaling up 
the 'hardware', i.e., achieving more qubits, gates, and longer coherence times. Much progress 
is being made on developing the necessary technology for a large-scale photonic quantum 
compute itM. The influence of noise is perhaps the most serious consideration!^ and must 
be overcome using error- correct ion and fault-tolerant constructiond^^^. We note that an 
alternative promising path to efficient quantum simulators is to exploit controllable quantum 
systems that can be used to directly implement model Hamiltonians, thereby avoiding the 
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aforementioned resource intensive approximation techniques and error correctiorP^. 
I. METHODS SUMMARY 

We use the minimal ST0-3G basi^ for H2, consisting of one |ls)-type atomic orbital per 
atom. These two functions are combined to form the bonding (antibonding) molecular 
orbital^. The 4 corresponding single-electron molecular spin-orbitals are combined an- 
tisymmetrically to form the six two-electron configurations ($1— >$6) that form the basis 
for our simulation. Due to symmetry, the Hamiltonian is block-diagonal in this basis, with 
blocks acting on each of the four subspaces spanned by {|$i), l^e)}) {|'^2)},{|'^'3)5 1'^4)}; and 
{|$5)} (See Methods, Section A). Therefore, finding the eigenvalues of the two 2x2 sub- 
matrices in the Hamiltonian — H^^'^') and H^^'^'' — amounts to performing the FCI. Estimating 
the eigenvalues of 2x2 matrices is the simplest problem for the IPEA. 

We employ a propagator time step of t=l h/E\, (the hartree, £'h~27.21 eV, is the atomic 
unit of energy), chosen so that 0<Et/27ch<l. For our proof-of-principle demonstration, 
all necessary molecular integrals are evaluated classically (Methods, Section C) using the 
Hartree-Fock procedurd^Sl. We use these integrals to calculate the matrix elements of H and 
U, then directly decompose each U^P''^^ operator into a logic gate network. We decompose 
the U^P'i^ operators into a global phase and a series of rotations of the one-qubit Hilbert 
spac^^ 

U = e'"Ry{P)R,{j)Ry{-P), (2) 

where a, (3, and 7, are real angles. is achieved by replacing angles a and 7 with ja and 
j7 (while P remains unchanged). Our decomposition of the controlled-t/-^ is shown in Fig.[l|D. 
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II. METHODS 



A. Minimal basis and symmetries in the electronic Hamiltonian of the hydrogen 
molecule 

The two |ls)-type atomic orbitals are combined to form the bonding and antibonding molec- 
ular orbital^, \g) and \u). Taking into account electron spin, the single-electron molecular 
spin-orbitals are denoted, \g]), \g[), \u]) and \u[), where 1 1) and | [) are the electron spin 
eigenstates. These are combined antisymmetrically to form the six two-electron configu- 
rations that form the basis for our simulation: \^i) = \g ] ,g [ \ = i\g] , 9 — g])) / ^/i!^)) 
\^^) = \g]^u] I, |$3) = |(7t,ni I, \^^) = \glu] I, \^^) = \glu[ \ and \^^) = \u],u[ \. Due 
to symmetry, the Hamiltonian is block-diagonal in this basis, with blocks acting on each of 
the four subspaces spanned by l^e)}, {|'^'2)}, {|*&3), |'^'4)}, and {l^s)}. Most of the 

elements of this basis are not mixed by the Hamiltonian. In particular, and l^g) niix 
only with each other because they have g symmetry while the rest have u symmetry. Of the 
remaining states only |$3) and |$4) mix because they have the same total z-projection of the 
spin, m5=0. [$2) and \(^^) have, respectively, ms=l and ms= — 1. Therefore, the Hamilto- 
nian is block- diagonal within four subspaces spanned by l^e)}; {|'^2)}, {l^s); 1^4)}) 
and {|$5)}. There are no approximations involved here, and finding the eigenvalues of the 
two 2x2 sub-matrices in the Hamiltonian (H^^'^^ and H^^'^^) amounts to performing an 
exact calculation (FCI) in the minimal basis. One should also note that it follows from 
the requirement that the wave functions are spin eigenstates, that the eigenstates of the 
subspace {|$3), |$4)} will be (|$3)±|$4))/v^. Additionally, there will be a three-fold de- 
generacy of the triplet state with angular momentum S=l. That is, the states |$2), l^s), 
and (|$3) + |$4))/v^ are degenerate. 

B. Details of computational methods 

Restricted Hartree-Fock calculations were carried out on a classical computer using the STO- 
3G basi^. The software used was the PyQuante quantum chemistry package version 1.6. 
The molecular integrals from the Hartree-Fock procedure are used to evaluate the matrix 
elements of the Hamiltonians H^^'^^ and H^^'^\ described in the main text. 

C. Classical error correction technique 

When running the IPEA, the probability of correctly identifying any individual bit with 
a single sample (^=1) is reduced from unity by both theoretical (inexact phase expansion 
to K bits) and experimental factors (such as imperfect gates). However, as long as it 
remains above 0.5, repeated sampling and a majority vote will reduce the probability of 
error exponentially with n, in accordance with the Chernoff bouncP^. This technique allows 
for a significant increase in success probability, at the expense of repeating the experiment 
a fixed number of times. We note that this simple classical error correction technique can 
only play a small role when it comes to dealing with errors in large-scale implementations. 
Here, the numerous errors in very large quantum logic circuits will make achieving a bit 
success probability over 0.5 a significant challenge, that must be met with quantum error 
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correction technique d^ * ^^ l 



D. Count rates 

We operate with a low-brightness optical source (spontaneous parametric downconversion 
pumping power ^50 mW) to reduce the effects of unwanted multi-photon-pair emissions 
(which cannot be distinguished by our non-photon-number-resolving detectors and intro- 
duce error into the circuit operation). This yields about 15 coincident detection events 
per second at the output of our optical circuit. Therefore each iteration can be repeated 
15 times a second. Reconfiguring the circuit for different iterations takes approximately 
7 seconds, largely due to the finite time required to rotate standard computer controlled 
waveplate mounts. Therefore, obtaining a 20-bit estimation of a phase takes about 3 min- 
utes, when using n=31 samples to determine the logical state of each bit (as was employed 
to achieve the results shown in Fig. |2]). Note that approximately 95% of this time is spent 
rotating waveplates. In future implementations, this time could be reduced significantly 
using integrated-photonics, e.g. qubit manipulation using an electrooptically-controUed 
waveguide Mach-Zehnder interferometei'^. 
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III. APPENDIX 



A. Efficient simulation of arbitrary molecular time-evolution operators 

A fundamental challenge for the quantum simulation of large molecules is the accurate de- 
composition of the system's time evolution operator, If. In our experimental demonstration, 
we exploit the small size and inherent symmetries of the hydrogen molecule Hamiltonian 
to implement U exactly, using only a small number of gates. As the system size grows 
such a direct decomposition will no longer be practical. However, an efficient first-principles 
simulation of the propagator is possible for larger chemical system^^ 

The key steps of an efficient approach are: (1) expressing the chemical Hamiltonian in 
second quantized form, (2) expressing each term in the Hamiltonian in a spin 1/2 rep- 
resentation via the Jordan- Wigner transformatiorP^, (3) decomposing the overall unitary 
propagator, via a Trotter-Suzuki expansiorpESI^ into a product of the evolution operators for 
non-commuting Hamiltonian terms, and (4) efficiently simulating the evolution of each term 
by designing and implementing the corresponding quantum circuit. We note that the first 
two steps generate a Hamiltonian that can be easily mapped to the state space of qubits. 
The last steps are part of the quantum algorithm for simulating the time-evolution operator, 
If, generated by this Hamiltonian. Details of each step are provided as follows: 



Step 1. Second- quantized Hamiltonian 

The general second-quantized chemical Hamiltonian has 0(A^^) terms, where N is the num- 
ber of single-electron basis functions (i.e. spin-orbitals) used to describe the system?^. The 
Hamiltonian can be explicitly written as: 

H = ^ ^ hpgdpa,q -\- — ^ ^ hpgrsO,p(i^Clr(ls, (SI) 
p,q p,q,r,s 

where the annihilation and creation operators {aj and respectively) obey the fermionic 
anti-commutation relations: [dj,a^]+ = 6ij and [ai,dj]+ = 0, and the indices p, q, r, and 
s run over all single-electron basis functions. The integrals hpg and hpgrs are evaluated 
during a preliminary Hartree-Fock procedur^^Hl g^^^j g^re defined as 

and 

U [a A X;(xi)Xg(x2)Xr(x2)x.(xi) 

hpgrs = / dxidx2-^^ 

J 

where Xg(x) are a selected single-particle basis. Here is the Laplacian with respect to 
the electron spatial coordinates, while r^x and ri2 are the distances between the a^^ nucleus 
and the electron and the distance between electrons 1 and 2, respectively. 
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Expressing the Hamiltonian in second-quantized form allows straightforward mapping of 
the state space to qubits. The logical states of each qubit are identified with the fermionic 
occupancy of a single-electron spin-orbital (i.e. |0) = occupied, |1) = unoccupied). There- 
fore, simulating a system with a total of single-electron spin-orbitals (e.g., N = Xk for 
a molecule with A atoms each with k, spin-orbitals) requires only qubits. Note that the 
A^-qubit Hilbert space allows for any number of electrons (up to A^), hence the scaling is in- 
dependent of the number of electrons present in the system. In practical Gaussian basis-set 
calculations, the number of spin-orbitals per atom is usually constant for a given row of the 
periodic tabld^. The use of a double-zeta basis setP^ would require employing ^ 30 logical 
qubits per simulated atom. For example, 1800 logical qubits would be required to store the 
wave function of the fuUerene (Cgo) molecule. 



Step 2. Jordan- Wigner transformation of the fermionic operators to spin variables 



Starting with the second-quantized Hamiltonian from (SI), the Jordan- Wigner transfor- 
mation is used to map fermionic creation and annihilation operators into a representation 
using the Pauli spin matrices as a basi^. This allows for a convenient implementation 
on a quantum computei *^^ * ^^ l The representation is achieved via the following invertible 



transformations, which are applied to each term in (SI): 



a, ^ l^^'-i ® a+ ® {^a'f^-i (S2a) 
a+ ^ 1^^-^ ®a-® (a^)®^-^' , (S2b) 

where a+ = (a^ + ia^)/2 = |0)(1| and a' = -iay)/2 = |1)(0|. The operators 
achieve the desired mapping of occupied (unoccupied) states to the computational basis 
[i.e., |1) (|0))] while other terms serve to maintain the required anti-symmetrization of the 
wavefunction in the spin (qubit) representation. 



Step 3. Exponentiation of the Hamiltonian 



As the system size represented by the chemical Hamiltonian (SI) grows, a direct de- 
composition of the time-evolution operator. If, into a sequence of logic gates will no longer 
be practical as the best methods scale exponentially. However, the Hamiltonian is a sum 
of one and two-electron terms whose time-evolution operators can each be implemented 
efficiently — e.g. with a number of gates that does not scale with A^. However, generally 
the terms do not commute, thus simple reconstruction of If from direct products of the 
individual operators is not possible. Trotter-Suzuki relations can be used to approximate 



the full unitary propagator from the individual evolution of non-commuting operators^'^. 
For a Hamiltonian H = J2i=i first-order Trotter-Suzuki decomposition is expressed 

as 

_t_ 

JJ(t) = e"*^* = (^e-''^''^'e-''^^'^' ■ ■ ■ e"^'^^'^*) * + 0{dt^). (S3) 

The value T„ = t/rft is called the Trotter numbei^^. As the Trotter number tends to infin- 
ity, or equivalently dt ^ 0, the approximation becomes exact. In practice, a compromise 
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between computational effort and accuracy is employed. In numerical computations, suc- 
cessive calculations at different timesteps dt are often carried out, and an extrapolation of 
dt ^ gives an estimate of the exact answer. A similar approach can be used for quantum 
simulation. 

We note that, unlike our small-scale experiment, the powers of the system evolution op- 
erator, U\ required for the IPEA cannot be achieved by simply changing parameters in the 
gate decomposition for U. In general U"^ will take twice as many gates as U. Intuitively, 
the system dynamics must be propagated for twice as long leading to twice as many ma- 
nipulations of the quantum simulator's natural dynamics. The increase in the number of 
gates required for extra bits will clearly amplify experimental errors, thereby limiting the 
obtainable precision. Note that although the number of required gates increases exponen- 
tially with the number of bits, each additional bit itself provides an exponential increase in 
precision. 

As mentioned in the manuscript, quantum algorithms that circumvent the problems found 
from the Trotter expansion are a fertile area of research. In the current scheme, Hamiltonians 
that are diagonal in the computational basis, such as the classical Ising model do not require 
a Trotter expansion for their accurate simulatiorP. 



Step 4- Circuit representations of the unitary propagator 



Each exponentiated tensor product of Pauli matrices can then be implemented efficiently 
by employing a family of quantum circuits. In order to provide an accurate estimation of 
an upper bound of the number of gates required for the different kinds of second-quantized 
operators, we carried out analytical gate decompositions. The circuit networks obtained 

The networks shown realize the unitary operator U{dt) for a 
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are summarized in Fig. 
general molecular Hamiltonian. To realize a controlled unitary, c — U{dt), as required by 
the phase estimation algorithm, only the rotations Rz{0) must be converted to controUed- 
Rz{0) rotations. The number of gates required to simulate each term is linear in the number 
of intervening qubits due to the product of terms resulting from the Jordan- Wigner 



transformation of Eq. S2 Therefore, the scaling of the number of quantum gates required 
for simulating a general many-electron chemical Hamiltonian is 0{N^) without considering 
the influence of noisd^. Fault tolerant quantum simulatioi]^ requires the use of a finite 
set of gates and the conversion from the continuous set of gates to a discrete set can be 
accomplished with polylogarithmic overheacP^. The encoding of robust quantum states 
will also require several redundant qubits for each logical qubit needed^. A more detailed 
analysis of fault tolerance in the context of quantum simulation can be found in Ref.l^. 



Resource count for a simple example. 

In order to illustrate this algorithm, we performed numerical simulations for H2 in the same 
minimal basis (ST0-3G) employed in our experiment. Unlike our experimental mapping, 
the logical states of each register qubit are now identified with the fermionic occupancy of 
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the four single-electron spin-orbitals (i.e. |0) = occupied, |1) = unoccupied). Therefore, the 
calculation requires a total of five qubits taking into consideration the single control qubit 
required for the IPEA. If quantum error correction is needed, the number of qubits will 
increase according to the scheme usecP^. Fig. S2 shows the error in the ground state energy 
as a function of the Trotter step. The ground state energies of the approximate unitary 
propagators were obtained via direct diagonalization on a classical computer. A precision of 
±10~'^Eh is achieved at a Trotter number of 6, which corresponds to 522 gates. Note that this 
gate count is to construct and includes both one- and two-qubit operations. This estimate 
does not take into consideration error correction for the qubits and it uses a continuous set 
of gates. In the path to large scale implementations, both will be serious considerations and 
will increase the complexity of the algorithm and the number of qubits necessarj^^*^. The 
unitary matrix must be raised to various powers to perform phase estimation. If one desires 
to maintain a fixed accuracy of 13 bits, about 8.5 x 10^ gates must be used for the IPEA 
estimation procedure. Note that this can be achieved by repeating the 522 gates required 
for U many times. Note that this does not include the resources associated with preparing 
a system eigenstate. If one uses an adiabatic state preparation technique^' the resources are 
proportional to the gap between the ground state and the excited state along the path of 
adiabatic evolutiorP^. 

Although the estimates just given exceed the capabilities of current quantum computers, 
these resource requirements grow only polynomially with the size the of system. Conse- 
quently, for large enough chemical systems, quantum computers with around 100 qubits are 
predicted to outperform classical computational devices for the first-principles calculation 
of chemical propertie j^ * ^° l 



B. Additional experimental results 

We use the estimation of the ground state energy at the equilibrium bond length, 1.3886 Qq 
(where oq is the Bohr radius) to study the effect of varying a range of experimental pa- 
rameters on the IPEA success probability, where we define the latter as the probability of 



correctly obtaining the phase to an precision of 2"*". Fig. S3i shows results measured over 
a range of n, the number of samples used to determine each bit. 

The probability of correctly identifying any individual bit with a single sample {n = 1) 
is reduced from unity by both theoretical {6) and experimental factors (such as imperfect 
gates). However, as long as it remains above 0.5, repeated sampling and a majority vote will 
improve the probability of correct identification. The data show that this is achieved and the 
error probability decreases exponentially with n, in accordance with the Chernoff bouncP^. 
This technique allows for a significant increase in success probability, at the expense of 
repeating the experiment a fixed number of times. We note that this simple classical error 
correction technique can only play a small role when it comes to dealing with errors in 
large-scale implementations. Here, the numerous errors in very large quantum logic circuits 
will make achieving a bit success probability over 0.5 a significant challenge, that must be 
met with quantum error correction techniqued^HSl^ 

Fig. [S3b shows the algorithm success probability measured as a function of the num- 
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ber of extracted bits (phase precision). By employing n = 101 samples per bit we achieve 
near perfect algorithm success probability up to 47 bits (yielding an energy precision of 

10~^^i?h), where this limit is imposed only by the machine-level precision used for the 
classical preprocessing of the Hamiltonians. It is insightful to understand how achieving such 
high precision will become a far more significant challenge for large-scale implementations: 
due to the small-scale of our demonstration, we are able to implement each power of U^'^'^^ 
directly, by re-encoding the same number of gates. Therefore, the probability of error intro- 
duced by gate imperfections remains a constant for each bit (and, in our implementation, 
under 50%). This is the main algorithmic feature that allows the high precision obtained 
in this experiment. However, as expounded in the appendix (section A), this will not be 
possible for larger implementations. In general, If will not have the same form as U"'. For 
each additional digit of precision sought, the gate requirements of the algorithm are roughly 
doubled, thereby amplifying any gate error. 

Fig. [S3]: shows the algorithm success probability measured as a function of the fidelity F 
(see caption) between the encoded register state and the ground state. The results show that 
our implementation is robust for F > 0.5. Because the probability of correctly obtaining 
each bit in a single measurement {n = 1) is greater than 0.5 in this regime, multiple sampling 
(n > 1) enables the success probability to be amplified arbitrarily close to unity. This is a 
general feature that will hold for large-scale implementations. However, for F < 0.5, the 
measured success probabilities are very low. 

If the register state output after each iteration is used as the input of the next, then the 
problem with low eigenstate fidelities can be overcome as the measurement of the control 
qubit collapses the wave function. Any pure encoded register state can be written in the 
eigenstate basis as \G) = ^jOjlAj), where is the fidelity of \G) with eigenstate jAj). 
Successful measurement of the m^^ bit associated with |Aj) will cause the register wavefunc- 
tion to collapse into a state with a greater fidelity with |Aj) — those eigenstates with a low 
probability of returning the measured bit value will be diminished from the superposition. 
As more bits are successfully measured, the register state will rapidly collapse to |Aj). In 
this way, the algorithm will return all the bits associated with |Aj) with probability at least!^ 
|a;i|^(l — e). With current technology, correct operation of our optical circuit requires de- 
structive measurement of both the control and register qubits after each IPEA iteration. 
Therefore, in our experiment the register state must be re-prepared for each iteration. 



1. How we obtain IPEA success probabilities 

Denoting the first m binary bits of a phase as = 0.(f)i(j)2---4'm, there is, in general, a 
remainder < 6 < 1, such that = + 52^™. To achieve an accuracy of ±2""^ the IPEA 
success probability is the sum of the probabilities for obtaining and + 2"*". This can 
be estimated experimentally, for a given phase, by simply repeating the algorithm a large 
number of times and dividing the number of acceptable results by the total. An estimate 
with an error less than 10% would require over 100 algorithm repetitions. We calculate the 
result shown in Fig. [SSb in this way. However, using this technique to obtain Fig. [SSb-c, and 



Fig. S3 (described below), would take a long time — the 20 points shown in each would require 
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more than 100 hours of waveplate rotation time alone. Instead, to obtain these results we 
force the appropriate feedforward trajectory {R{uJk)) for each accepted phase value and use 
n = 301 samples to estimate the 0/1 probabilities for each bit. Using the standard binomial 
cumulative distribution function it is then possible to calculate the majority vote success 
probability for each bit of each accepted value for a given n (1 and 101 in the figures). The 
probability for obtaining an accepted phase value is then the product of the majority vote 
success probabilities for each bit, and the total algorithm success probability is the sum of the 
probabilities for obtaining each accepted phase. The error bars represent a 68% confidence 
interval and are obtained from a direct Monte-Carlo simulation of the above process. 

Note that forcing the correct feedforward in this way, and taking many samples to estimate 
the 0/1 probabilities for each bit, simply allows us to accurately estimate the probability 
that the algorithm will return the correct phase by itself - i.e. without forcing the correct 
feedforward. 



2. Experimental model 



A simple computational model of our experiment produced the lines shown in Fig. [S3j This 
model allows for two experimental imperfections, which are described below, but otherwise 
assumes perfect optic element operation. The model consists of a series of operators, rep- 
resenting optical elements and noise sources, acting on a vector space representing both 
photonic polarisation and longitudinal spatial modd^. Firstly the model allows for photon 
distinguishability, quantified by an imperfect relative non-classical interference visibility of 
0.93 (ideal 1), which reduces the quality of our two-qubit logic gate. Secondly the model 
allows for phase damping of the control qubit, described by the operation element^^ 







v^r^ 



and 







(S4) 



Our model employs 7 = 0.06 (ideal 0), which corresponds to ~ 3% dephasing. These 
experimental imperfections are attributed to a combination of residual higher-order photon 
pair emissions from our optical source and circuit alignment drift during long measurement 
sets. 
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FIG. 1: Algorithm and experimental implementation, (a) IPEAEHH at iteration k. To 
produce a i^-bit approximation to (f) tlie algoritlim is iterated K times. Eacli iteration obtains one 
bit of (j) {(pk)'- starting from the least significant {4>k), k is iterated backwards from ii' to 1. The 
angle uj^ depends on all previously measured bits, uJk= ~ 2vr6, where b, in the binary expansion, 
is b=0.0cl)k+i4'k+2---4'K and ujk=0- H is the standard Hadamard gat^^. (b) Our gate network 
for a two-qubit controlled-C/-' gate, as discussed in the Methods Summary, (c) Two-qubit optical 
implementation of (a). Photon pairs are generated by spontaneous parametric down-conversion 
(SPDC), coupled into single-mode optical fiber and launched into free space optical modes C 
(control) and R (register). Transmission through a polarizing beamsplitter (PBS) prepares a 
photonic polarization qubit in the logical state |0), the horizontal polarization. The combination 
of a PBS with half (A/2) and quarter (A/4) waveplates allows the preparation (or analysis) of an 
arbitrary one-qubit pure state. The optical controlled- i?^ gate, shown in the dashed box, is realized 
using conditional transformations via spatial degrees of freedom as described by LanyoiP^ et al. 
Coincident detection events (3.1 ns window) between single photon counting modules (SPCM's) 
Dl and D3 (D2 and D3) herald a successful run of the circuit and result (1) for cpk- Waveplates 
are labelled with their corresponding operations. 
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FIG. 2: Quantum algorithm results: H2 potential energy curves in a minimal basis. Each 
point is calculated using a 20-bit IPEA and employing n=31 samples per bit (repetitions of each 
iteration). Every case was successful, achieving the target precision of ib(2~^^x27r) -E'h~10~^ E]^. 
Curve G (E3) is the low (high) eigenvalue of H^^'^\ Curve El is a triply degenerate spin-triplet 
state, corresponding to the lower eigenvalue of H^^'^^ as well as the eigenvalues H^'^^ and H^^\ 
Curve E2 is the higher (singlet) eigenvalue of H^^'^K Measured phases are converted to energies E 
via E=2TT(p+l/r, where the last term accounts for the proton-proton Coulomb energy at atomic 
separation r, and reported relative to the ground state energy of two hydrogen atoms at infinite 
separation. Inset a): Curve G rescaled to highlight the bound state. Inset b): Example of 
raw data for the ground state energy obtained at the equilibrium bond length, 1.3886 a.u.. The 
measured binary phase is (/>=0. 01001011101011100000 which is equal to the exact value, in our 
minimal basis, to a binary precision of ±2"^^. Note that the exact value has a remainder of (5~0.5 
after a 20 bit expansion, hence the low contrast in the measured 20*'^ bit. 
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FIG. SI: The quantum circuits corresponding to evolution of the listed Hermitian 
second-quantized operators. Here p, q, r, and s are orbital indices corresponding to qubits 
such that the population of |1) determines the occupancy of the orbitals. It is assumed that 
the orbital indices satisfy p > q > r > s. These circuits were found by performing the Jordan- 
Wigner transformation given in ( |S2b ) and (S2a) and then propagating the obtained Pauli spin 
variableJ^. In each circuit, 6 = 0{h) where h is the integral preceding the operator. Gate T{9) 
is defined by T|0) = |0) and r|l) = exp(— G is the global phase gate given by exp(— 
and the change-of-basis gate Y is defined as 7r/2). Gate H refers to the Hadamard gate. For 
the number-excitation operator, both M = Y and \\A = H must be implemented in succession. 
Similarly, for the double excitation operator each of the 8 quadruplets must be implemented in 
succession. The global phase gate must be included due to the phase-estimation procedure. Phase 
estimation requires controlled versions of these operators which can be accomplished by changing 
all gates with ^-dependence into controlled gates. 
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FIG. S2: Trotter error analysis and resource count for hydrogen molecule using a 
scalable quantum simulation algorithm, (a) Plot of ground state energy of hydrogen molecule 
as a function of the length of the time step. As the time step length decreases, the accuracy of 
the approximation increases in accordance with eqn. (S3). The total time of propagation, t, was 
unity and this time was split into time steps, dt. The circles are at integer values of the Trotter 
number, Tn = t/dt. Green horizontal lines indicate the bounds for ±W~^Eh precision, (b) Gates 
for a single construction of the approximate unitary as a function of time step. As the time step 
decreases, more gates must be used to construct the propagator. The triangles indicate integer 
values of the Trotter number and the green vertical line corresponds to the same threshold from 
graph a. Perfect gate operations are assumed. 
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FIG. S3: IPEA success probability measured over a range of parameters. Probabilities 
for obtaining tlie ground state energy, at tlie equilibrium bond length 1.3886 oq, as a function 
of: (a) the number of times each bit is sampled (n); (b) the number of extracted bits (m); (c) 
the fidelity between the encoded register state and the ground state (F). The standard fidelitjff^ 
between a measured mixed p and ideal pure state is -F=(^'|p|^'). (a) & (b) employ a ground 
state fidelity of F ~ 1. (a) & (c) employ a 20-bit IPEA. All lines are calculated using a model 
that allows for experimental imperfections. This model, as well as the technique used to calculate 
success probabilities and error bars, are detailed in the appendix (section B). 



20 



